TFM Caracterización inmunonómica del microambiente tumoral en cáncer colorrectal

---

1 Obtención y preparación de los datos

Cargamos los arhivos contenidos en “infoClin_all_dataset.RData” y “CRC_stageII.RData” con los contenidos clínicos y de CRC.

# Cargamos el archivo infoClin_all_dataset.RData 16/03/2024
dataClin <- load("data/infoClin_all_dataset.RData")

# Cargamos el archivo CRC_stageII.RData 19/03/2024 dataCRC <-
# load('data/CRC_stageII.RData')

# Cargamos el archivo exp_mat.RData 02/04/2024
dataCRC <- load("data/exp_mat.RData")

GSE Files

• GSE39582 (Marisa et al. 2013) 421 vs 585
• GSE14333 (Jorissen et al. 2009) 185 vs 290
• GSE13294 (Jorissen et al. 2008) 121 vs 155
• GSE17536 (Smith et al. 2010) 111 vs 177

Platform GPL570: [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array

Fusionamos los datos en un único dataframe al que le añadimos la columna grupo con los datos stage y msi.

Factorizamos las variables categóricas y ordenamos por ID tanto los datos clínicos como los de la matriz de expresión.

Creamos un DataFrame para asignar cada ID a su estudio clínico. A cada estudio le asignamos un color para poder identificarlo fácilmente en los estudios gráficos.

# 11/04/2023 Función para preparar el data frame con el estudio por ID. También
# asignaremos un color a cada estudio. Utilizamos el paquete 'dplyr'
prepare_df <- function(df, study_name, color) {
    df %>%
        select(ID) %>%
        mutate(study = study_name, color = color)
}

# Lista de estudios con sus colores asociados
studies_info <- list(clx = "red", GSE39582 = "blue", GSE14333 = "green", GSE13294 = "orange",
    GSE17536 = "yellow", TCGA = "violet")

# Creamos una lista de de los estudios utilizando la función prepare_df
df_studies_list <- lapply(names(studies_info), function(x) {
    prepare_df(get(paste0(x, "_infoClin")), x, studies_info[[x]])
})

# Juntamos los df's de la lista se studies en uno solo y lo ordenamos por 'ID'
df_studies <- bind_rows(df_studies_list) %>%
    arrange(ID)

# Fusionamos y ordenamos todos los datos clínicos generando 'df_dataClin'
df_dataClin <- bind_rows(lapply(names(studies_info), function(x) get(paste0(x, "_infoClin")))) %>%
    arrange(ID)

# Cramos una nueva columna 'grupo' a los 'df_dataClin' con los datos de msi y
# stage
df_dataClin$grupo <- with(df_dataClin, paste(stage, msi_imputed, sep = "_"))

# Cramos una nueva columna 'cms_stage' a los 'df_dataClin' con los datos de cms
# y stage 19/04/2024
df_dataClin$cms_stage <- with(df_dataClin, ifelse(is.na(cms), NA, paste(stage, cms,
    sep = "_")))

# Convertimos a factores
df_dataClin <- df_dataClin %>%
    mutate(across(c(msi_imputed, cms, stage, grupo, cms_stage), factor))

# Añadimos la columna 'study' a df_dataClin. Al final no añadimos 'color'
df_dataClin <- left_join(df_dataClin, select(df_studies, ID, study), by = "ID")

Pasamos la matriz de expresión a data_frame y la ordenamos por “ID”

# Pasamos la matriz de expresión a data frame 26/03/2024
df_exmat <- as.data.frame(ex)

# Ordenamos la matriz de expresión por nombre de columnas 26/03/2024
df_exmat <- df_exmat[, order(names(df_exmat))]

# Transponemos la matriz de expresión: 22/04/2025
df_exmat_t <- as.data.frame(t(df_exmat))

# Creamos un df con los datos clínicos y los de expresión génica (transpuesta)
df_dataClin_exmat_t <- cbind(df_dataClin, df_exmat_t)

2 Exploración de los datos

Resumen exploratorio de los datos:

# Hacemos un Summary con el paquete skimr 18/03/2024
skim(df_dataClin)

Data summary

Name	df_dataClin
Number of rows	1079
Number of columns	7
_______________________
Column type frequency:
character	2
factor	5
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
ID	0	1	7	12	0	1079	0
study	0	1	3	8	0	6	0

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
msi_imputed	0	1.00	FALSE	2	MSS: 840, MSI: 239
cms	158	0.85	FALSE	4	CMS: 329, CMS: 262, CMS: 194, CMS: 136
stage	0	1.00	FALSE	2	II: 684, III: 395
grupo	0	1.00	FALSE	4	II_: 511, III: 329, II_: 173, III: 66
cms_stage	158	0.85	FALSE	8	II_: 210, II_: 151, II_: 124, III: 119

Estadísticas y cálculos para la generación de los gráficos de datos clínicos:

# Estadísticas desglosadas 18/03/2024 Frecuencias
msi_frecuencias <- table(df_dataClin$msi_imputed)
cms_frecuencias <- table(df_dataClin$cms)
stage_frecuencias <- table(df_dataClin$stage)

# Valores NA
na_msi <- sum(is.na(df_dataClin$msi_imputed))
na_cms <- sum(is.na(df_dataClin$cms))
na_stage <- sum(is.na(df_dataClin$stage))

# Contamos cuantos MSI y MSS tenemos por cada estadio, II y III 19/03/2024
sII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "II")
sII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "II")
sIII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "III")
sIII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "III")

# Crear un dataframe con estos valores
stage_MSI_MSS_df <- data.frame(Stage_MSI_MSS = c("II_MSI", "II_MSS", "III_MSI", "III_MSS"),
    Freq = c(sII_MSI, sII_MSS, sIII_MSI, sIII_MSS))

# Contamos cuantos MSI y MSS tenemos por cada estadio, II y III 19/04/2024
sII_CMS1 <- sum(df_dataClin$cms_stage == "II_CMS1", na.rm = TRUE)
sII_CMS2 <- sum(df_dataClin$cms_stage == "II_CMS2", na.rm = TRUE)
sII_CMS3 <- sum(df_dataClin$cms_stage == "II_CMS3", na.rm = TRUE)
sII_CMS4 <- sum(df_dataClin$cms_stage == "II_CMS4", na.rm = TRUE)
sIII_CMS1 <- sum(df_dataClin$cms_stage == "III_CMS1", na.rm = TRUE)
sIII_CMS2 <- sum(df_dataClin$cms_stage == "III_CMS2", na.rm = TRUE)
sIII_CMS3 <- sum(df_dataClin$cms_stage == "III_CMS3", na.rm = TRUE)
sIII_CMS4 <- sum(df_dataClin$cms_stage == "III_CMS4", na.rm = TRUE)

# Crear un dataframe con estos valores
stage_CMS_df <- data.frame(Stage_CMS = c("II_CMS1", "II_CMS2", "II_CMS3", "II_CMS4",
    "III_CMS1", "III_CMS2", "III_CMS3", "III_CMS4"), Freq = c(sII_CMS1, sII_CMS2,
    sII_CMS3, sII_CMS4, sIII_CMS1, sIII_CMS2, sIII_CMS3, sIII_CMS4))

3 Análisis de los datos clínicos

3.1 ESTADIO

# Pasamos a DataFrame stage
df_stage_frecuencias <- as.data.frame(stage_frecuencias)
df_stage_frecuencias

  Var1 Freq
1   II  684
2  III  395

# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_stage_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_stage_frecuencias$Porcentaje <- (df_stage_frecuencias$Freq / total_frecuencias) * 100

# Diagrama de Barras de ESTADIO
p1 <- ggplot(df_stage_frecuencias, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
  geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
  labs(title = "STAGE", x = "", y = "Frecuencia") +
  theme_minimal()
p1

3.2 MSI/MSS

# Pasamos a DataFrame msi
df_msi_frecuencias <- as.data.frame(msi_frecuencias)
df_msi_frecuencias

  Var1 Freq
1  MSI  239
2  MSS  840

# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_msi_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_msi_frecuencias$Porcentaje <- (df_msi_frecuencias$Freq / total_frecuencias) * 100

# Diagrama de Barras de MSI/MSS
p2 <- ggplot(df_msi_frecuencias, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
  geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
  labs(title = "MSI/MSS", x = "", y = "Frecuencia") +
  theme_minimal()
p2

3.3 MSI/MSS & ESTADIO

# Stage vs MSI/MSS DataFrame
stage_MSI_MSS_df

  Stage_MSI_MSS Freq
1        II_MSI  173
2        II_MSS  511
3       III_MSI   66
4       III_MSS  329

# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(stage_MSI_MSS_df$Freq)
# Añadimos una columna de porcentajes al dataframe
stage_MSI_MSS_df$Porcentaje <- (stage_MSI_MSS_df$Freq / total_frecuencias) * 100

# Diagrama de Barras de Stage vs MSI/MSS
p3 <- ggplot(stage_MSI_MSS_df, aes(x = Stage_MSI_MSS, y = Freq)) +
  geom_bar(stat = "identity", fill = c("salmon", "seashell3",
                                       "turquoise3", "slategray")) +
  geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 2, size = 3.5) +
  labs(title = "STAGE & MSI/MSS", x = "", y = "Frecuencia") +
  theme_minimal()
 
p3

3.4 CMS

# Pasamos a DataFrame cms
df_cms_frecuencias <- as.data.frame(cms_frecuencias)
df_cms_frecuencias

  Var1 Freq
1 CMS1  194
2 CMS2  329
3 CMS3  136
4 CMS4  262

# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_cms_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_cms_frecuencias$Porcentaje <- (df_cms_frecuencias$Freq / total_frecuencias) * 100

# Diagrama de Barras de CMS
p4 <- ggplot(df_cms_frecuencias, aes(x = Var1, y = Freq, fill = Var1)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Freq), vjust = -0.3, size = 3.5) +  # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +  
  scale_fill_manual(values = c("salmon", "seashell3", "turquoise3", "slategray")) +
  labs(title = "CMS", x = "", y = "Frecuencia") +
  theme_minimal() +
  theme(legend.position = "none") # Omitir la leyenda

p4

3.5 CMS & ESTADIO

# Stage vs MSI/MSS DataFrame
stage_CMS_df

  Stage_CMS Freq
1   II_CMS1  124
2   II_CMS2  210
3   II_CMS3   96
4   II_CMS4  151
5  III_CMS1   70
6  III_CMS2  119
7  III_CMS3   40
8  III_CMS4  111

# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(stage_CMS_df$Freq)
# Añadimos una columna de porcentajes al dataframe
stage_CMS_df$Porcentaje <- (stage_CMS_df$Freq / total_frecuencias) * 100

# Diagrama de Barras de Stage vs MSI/MSS
p5 <- ggplot(stage_CMS_df, aes(x = Stage_CMS, y = Freq)) +
  geom_bar(stat = "identity", fill = c("salmon", "seashell3",
                                       "turquoise3", "slategray", 
                                       "salmon", "seashell3",
                                       "turquoise3", "slategray")) +
  geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 2, size = 3.5) +
  labs(title = "STAGE & CMS", x = "", y = "Frecuencia") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))# Inclinar las etiquetas del eje x
 
p5

3.6 TODOS

# Múltiple Layout con patchwork En un grid 3x2
p2 + p3 + p4 + p5 + p1 + plot_layout(ncol = 2, nrow = 3)

4 Análisis PCA

Para el Análisis de Componentes Principales (PCA), normalizamos y escalamos (estandarizamos) la matriz de expresión. En nuestro caso, los datos ya están normalizamos (en principio con \(log_2\)) por lo que tan solo los escalamos.
La estandarización se realiza restando la media de cada variable y dividiendo el resultado por la desviación estándar de esa variable.

\[ Z = \frac{x-\mu}{\sigma} \] Utilizamos la función scale para escalar. A continuación realizamos el PCA con la transpuesta de la matriz de expresión dado que queremos un análisis de las muestras en función de los genes y no al revés.
IMP: al hacer el escalado antes de la transposición estamos estandarizando por genes y no por muestras.

# Escalamos/normalizamos la matriz de expresión. También centra. 28/04/2024 Al
# hacer el escalado antes de la transposición estamos escalando por genes y no
# por muestras.
df_exmat_scaled <- scale(df_exmat)
# Hacemos el PCA. Hay que hacer la transpuesta de la matriz de expresión (las
# columnas representan genes y las filas representen muestras), así los
# componentes principales reflejan la variabilidad entre muestras. 28/04/2024
# Na haría falta center = TRUE ya que al escalar ya centralizamos.
pca_result <- prcomp(t(df_exmat_scaled), center = TRUE, scale. = FALSE)

RESUMEN PCA

# Obtenemos los valores propios (la varianza explicada por cada componente)
varianza <- pca_result$sdev^2

# Calculamos la proporción de varianza explicada por cada componente
prop_varianza <- varianza/sum(varianza)

# Calculamos la proporción acumulativa de varianza explicada
prop_acumulada <- cumsum(prop_varianza)

# Crear un dataframe con esta información para las primeras 10 PC
resumen_pca <- data.frame(SD = pca_result$sdev[1:10], Varianza = prop_varianza[1:10],
    VarAcumulada = prop_acumulada[1:10])

print(resumen_pca)

          SD   Varianza VarAcumulada
1  13.721545 0.08197744   0.08197744
2  11.124090 0.05387873   0.13585617
3  10.811433 0.05089263   0.18674880
4   9.600035 0.04012675   0.22687555
5   6.928531 0.02090117   0.24777672
6   6.605919 0.01900005   0.26677677
7   6.062676 0.01600358   0.28278035
8   5.642600 0.01386267   0.29664302
9   5.370890 0.01255975   0.30920276
10  5.178365 0.01167545   0.32087821

La baja varianza explicada por los primeros dos componentes, PC1 (8,2%) y PC2 (5.4%), sugiere que los datos tienen una distribución de varianza muy dispersa y que no hay unas pocas direcciones principales de variabilidad que dominen.
Esto es coherente con el hecho de que estamos tratando con datos biológicos, en concreto de CRC, que suelen ser de alta dimensionalidad con múltiples factores contribuyendo a la variabilidad general. Es decir, no hay patrones dominantes de variabilidad que puedan ser fácilmente capturados por los primeros componentes principales.

CARGAS DE GENES EN PCA

Para entender mejor qué genes están contribuyendo a la separación observada, podemos analizar las cargas de los componentes principales para ver qué genes tienen mayores pesos en los componentes que diferencian los grupos.

# Extracción de cargas de componentes principales
cargas_pca <- pca_result$rotation

# Mostramos las cargas de los primeros dos componentes
primer_componente <- cargas_pca[, 1]
segundo_componente <- cargas_pca[, 2]

PC1

# Extraemos los top 10 genes y sus cargas de PC1
top_genes_pc1 <- head(sort(primer_componente, decreasing = TRUE), 10)

# Convertimos a data frame para visualización de PC1
top_genes_pc1_df <- data.frame(Gen = names(top_genes_pc1), Carga = top_genes_pc1)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc1_df <- top_genes_pc1_df[order(top_genes_pc1_df$Carga, decreasing = TRUE),
    ]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc1_df$Gen <- factor(top_genes_pc1_df$Gen, levels = top_genes_pc1_df$Gen)

# Creamos el gráfico de barras
plot_pc1 <- ggplot(top_genes_pc1_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
    position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
    hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC1", x = "Gen",
    y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")

# Mostramos los top 10 genes de PC1
print(top_genes_pc1)

     SFRP2       GAS1    COL10A1      CCL18     CYP1B1       SPP1      GPNMB 
0.06248246 0.04743630 0.04494118 0.04386073 0.04165178 0.03902618 0.03896303 
     GREM1     CCDC80      POSTN 
0.03860923 0.03848662 0.03844061 

PC2

# Extraemos los top 10 genes y sus cargas de PC2
top_genes_pc2 <- head(sort(segundo_componente, decreasing = TRUE), 10)

# Convertimos a data frame para visualización de PC2
top_genes_pc2_df <- data.frame(Gen = names(top_genes_pc2), Carga = top_genes_pc2)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc2_df <- top_genes_pc2_df[order(top_genes_pc2_df$Carga, decreasing = TRUE),
    ]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc2_df$Gen <- factor(top_genes_pc2_df$Gen, levels = top_genes_pc2_df$Gen)

# Creamos el gráfico de barras
plot_pc2 <- ggplot(top_genes_pc2_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
    position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
    hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC2", x = "Gen",
    y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")

# Mostramos los top 10 genes de PC1
print(top_genes_pc2)

      REG4       ZIC2      REG1A      TRIM7       CTSE       TCN1       AGR3 
0.07056229 0.05049736 0.04890627 0.04815804 0.04715623 0.04230700 0.04170686 
    RAB27B    SDR16C5    PLA2G2A 
0.04125286 0.04028539 0.03911398 

# Mostramos gráficamente las cargas del top 10 tanto de PC1 como de PC2 en un
# grid 2x1 (patchwork)
plot_pc1 + plot_pc2 + plot_layout(ncol = 2, nrow = 1)

GRÁFICAS PCA

Utilizamos el paquete factoextra para realizar las gráficas de los 2 primeros componentes, en concreto cla función fviz_pca_ind, que además nos permite añadir las elipses de confianza.

4.1 Estudios de origen

Entre los diferentes grupos de estudio, las elipses de confianza están centradas y solapadas lo que sugiere que, en el espacio de los primeros dos componentes principales, no hay una separación clara entre los grupos.

No es descartable la necesidad de examinar componentes adicionales o utilizar otras técnicas de análisis como métodos de Clustering.

# Utilizamos "factoextra" para hacer la gráfica de la PCA
# Hacemos la gráfica de los dos primeros componentes
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$study, # Colores de los puntos basados en los estudios
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "Study",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA Estudios")

4.2 MSI/MSS

Cuando representamos la PCA por inestabilidad de microsatélites vemos que las dos elipses, las de MSI y las de MSS, están claramente separadas lo que sugiere que ambos grupos tienen perfiles de expresión genética distintos en las dimensiones capturadas por los componentes principales (CP).
Como ya sabemos, estas diferencias son biológicamente relevantes; los componentes que muestran esta separación probablemente están capturando las variaciones genéticas ya conocidas de antemano.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$msi_imputed, # Colores de los puntos basados en MSI/MSS
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "MSI/MSS",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA MSI/MSS")

4.3 ESTADIO II/III

Entre los dos grupos de estadio, II y III, las elipses de confianza están más o menos centradas y algo solapadas lo que sugiere que, no hay una separación clara entre los grupos. Algo un tanto previsible ya que el estadio del cáncer depende del tiempo de progreso de la enfermedad.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$stage, # Colores de los puntos basados en Estadio
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "Estadio",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA Estadio II/III")

4.4 MSI/MSS vs ESTADIO II/III

Se podría decir que esta gráfica és una combinación de las dos anteriores; fijándonos en las elipses de confianza, se observan dos grupos separados claramente de otros dos, debido principalmente a la contribución de la instabilidad de microsatélites y no tanto a la del estadio del CRC.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$grupo, # Colores de los puntos basados en Grupos
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "Grupo",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA Grupos")

4.5 CMS

Al igual que con la inestabilidad de microsatélites, cuando representamos la PCA por clasificación CMS vemos que las cuatro elipses están claramente separadas; los 4 grupos tienen perfiles de expresión genética distintos en las dimensiones capturadas por los CP.
Estas diferencias son biológicamente relevantes dado que los componentes que muestran esta separación están capturando variaciones genéticas ya conocidas.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$cms, # Colores de los puntos basados en CMS
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "CMS",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA por CMS")

4.6 CMS-ESTADIO

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$cms_stage, # Colores de los puntos basados en CMS
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "cms_stage",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA por CMS-Estadio")

5 Estimación de la infiltración celular

Analizamos los datos viendo la infiltración celular.

Utilizamos el paquete MCPcounter (Becht et al. 2016) que estima la abundancia poblacional de poblaciones de células inmunes (8 tipos) y estromales (2 tipos) que se infiltran en tejidos mediante la expresión genética.

En concreto, utilizatemos MCPcounter para estimar la infiltración celular del estroma y el sistema inmunitario en el microambiente tumoral (TME).

# Cargamos la librería 'MCPcounter'
require(MCPcounter)

# Con todos los datos de la matriz de expresión 26/03/2024
exmat_MCP <- MCPcounter.estimate(expression = df_exmat, featuresType = "HUGO_symbols")
# Transponemos
exmat_MCP <- t_df(exmat_MCP)
# Añadimos los datos clínicos a los datos de la matriz generada con MCPcounter
# 26/03/2024
exmat_MCP_Clin <- cbind(df_dataClin, exmat_MCP)

# Guardamos 'df_dataClin', 'df_exmat' y 'exmat_MCP_Clin' en un archivo .Rdata.
# Lo hago para poder utilizar estos datos ya calculados en otros .rmd
# 21/04/2024 save(df_dataClin, df_exmat, exmat_MCP_Clin, file =
# 'tfm_SKF.RData')

A modo de recordatorio didáctico de Inferencia Estadística:

• El nivel de significación \(\alpha\) de un contraste es el error máximo de tipo I que estamos dispuestos a asumir.

	\(H_0\) es verdadera	\(H_0\) es falsa
Rechazar \(H_0\)	Error Tipo I (α)	Decisión Correcta
No rechazar \(H_0\)	Decisión Correcta	Error Tipo II (β)

• El p-valor es la probabilidad del resultado del estadístico de contraste observado o de uno más alejado cuando la hipótesis nula es cierta.

• El p-valor asociado a una observación del estadístico de contraste es el menor nivel de significación que nos permite rechazar la hipótesis nula.

• Si el p-valor es inferior al nivel de significación \(\alpha\), rechazaremos la hipótesis nula, en caso contrario la aceptaremos.

En nuestro caso utilizamos un nivel de significación \(\alpha = 0.05\)

5.1 Infiltración MSI/MSS vs stage

5.1.1 MSI/MSS vs stage de todos los tipos celulares

Para hacer una comparativa entre los 4 grupos, utilizaremos el test no paramétrico de Kruskal-Wallis. Se suele utilizar para determinar si hay diferencias significativas entre tres o más grupos independientes cuando los datos no cumplen los supuestos necesarios para realizar un ANOVA: normalidad y homocedasticidad (homogeneidad de varianzas).

# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024

# tests de Kruskal-Wallis para NK
kw_NK_result <- kruskal.test(exmat_MCP_Clin$`NK cells` ~ grupo, data = exmat_MCP_Clin)

# Para el resto de tipos celulares utilizamos una lista con los nombres
tipos_celulares <- c("NK cells", "T cells", "CD8 T cells", "Cytotoxic lymphocytes",
    "B lineage", "Monocytic lineage", "Myeloid dendritic cells", "Neutrophils", "Endothelial cells",
    "Fibroblasts")

# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw <- lapply(tipos_celulares, function(cell_type) {
    kruskal.test(reformulate("grupo", response = cell_type), data = exmat_MCP_Clin)
})

# Añadimos los nombres de los tipos celulares
names(resultados_kw) <- tipos_celulares

# Extraemos los valores de p de todos los tests realizados
p_valores <- sapply(resultados_kw, function(x) x$p.value)

La hipótesis nula del test de Kruskal-Wallis \(H_0\) establece que todas las poblaciones (grupos) tienen distribuciones idénticas, es decir, las medianas de todos los grupos son iguales.

\[ H_0: \forall i, j \, (i \neq j) \Rightarrow \mu_i = \mu_j \]

La hipótesis alternativa \(H_1\) establece que al menos uno de los grupos tiene una distribución diferente en cuanto a la mediana, comparada con las otras poblaciones sin especificar cuál de los grupos es diferente.

\[ H_1: \exists i \neq j \mid \mu_i \neq \mu_j \]

# Generamos todos los boxplot 27/03/2024
grupo <- exmat_MCP_Clin$grupo
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`NK cells`, grupo,
    "NK cells", p_valores["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`T cells`, grupo,
    "T cells", p_valores["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`CD8 T cells`, grupo,
    "CD8 T cells", p_valores["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
    grupo, "Cytotoxic lymphocytes", p_valores["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`B lineage`, grupo,
    "B lineage", p_valores["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Monocytic lineage`,
    grupo, "Monocytic lineage", p_valores["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Myeloid dendritic cells`,
    grupo, "Myeloid dendritic cells", p_valores["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Neutrophils, grupo,
    "Neutrophils", p_valores["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Endothelial cells`,
    grupo, "Endothelial cells", p_valores["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Fibroblasts, grupo,
    "Fibroblasts", p_valores["Fibroblasts"])

# En un grid 3x2 (patchwork)
bp1 + bp2 + bp3 + plot_layout(ncol = 3, nrow = 1)

bp4 + bp5 + bp6 + plot_layout(ncol = 3, nrow = 1)

bp7 + bp8 + bp9 + plot_layout(ncol = 3, nrow = 1)

bp10 + plot_layout(ncol = 3, nrow = 1)

# Todos juntos para el informe del TFM 11/06/2024
bp1 + bp2 + bp3 + bp4 + bp5 + bp6 + bp7 + bp8 + bp9 + bp10 + plot_layout(ncol = 4,
    nrow = 3)

5.1.2 Infiltración NKs

NORMALIDAD Y HOMOCEDASTICIDAD

Primero compromabos la normalidad de los grupos para las NK cells.

Tabla con los resultados de los tests de Shapiro para normalidad.
Si p-value < 0.05 rechazamos la hipótesis nula de igualdad de estadísticos.

Shapiro test	MSI	MSS	II	III
p-valor	3.339e-05	1.058e-10	1.371e-10	1.089e-04

Dado que no hay normalidad no hace falta comprobar homocedasticidad.

Para realizar una comparativa entre dos grupos en los que asumimos no normalidad , en nuestro caso MSI vs MSS y estadio II vs estadio III, utilizamos el test no paramétrico de Mann-Whitney U también conocido como el test de Wilcoxon en el que se asume que los datos no necesitan ser normalmente distribuidos aunque deben poder ordenarse.

# Para dos grupos utilizamos el test no paramético Mann-Whitney U (asumimos que
# no tenemos normalidad) 16/04/2024

# Realizamos el test de Mann-Whitney U para los grupos MSI y MSS
mwU_NK_msi_test <- wilcox.test(grupo_msi, grupo_mss, alternative = "two.sided")

# Realizamos el test de Mann-Whitney U para los estadios II y III
mwU_NK_stage_test <- wilcox.test(grupo_sII, grupo_sIII, alternative = "two.sided")

La hipótesis nula \(H_0\) para el test de Mann-Whitney U establece que no hay diferencia en las medianas (o, más generalmente, en las distribuciones) de los dos grupos: \[ H_0: \mu_{grupo1} = \mu_{grupo2} \] La hipótesis alternativa \(H_1\) establece que hay una diferencia entre las distribuciones de los dos grupos. En nuestro caso la planteamos Bilateral (Two-sided):

\[ H_1: \mu_{grupo1} \neq \mu_{grupo2} \]

# Generamos todos los boxplot referentes a NKs 27/03/2024 Diagrama de Barras de
# NKs MSI vs MSS
bp11 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`NK cells`,
    exmat_MCP_Clin$msi_imputed, "NK cells MSI vs MSS", mwU_NK_msi_test$p.value)
# Diagrama de Barras de NKs II vs III
bp12 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`NK cells`,
    exmat_MCP_Clin$stage, "NK cells II vs III", mwU_NK_stage_test$p.value)

# En un grid 2x2 (patchwork)
p4 + bp1 + p1 + bp11 + p3 + bp12 + plot_layout(ncol = 2, nrow = 3)

5.2 Infiltración resto de tipos celulares

5.2.1 T cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp13 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`T cells`,
    exmat_MCP_Clin$msi_imputed, "T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp14 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`T cells`,
    exmat_MCP_Clin$stage, "T cells II vs III")

# En un grid 1x3 (patchwork)
bp2 + bp13 + bp14 + plot_layout(ncol = 3, nrow = 1)

5.2.2 CD8 T cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp15 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`CD8 T cells`,
    exmat_MCP_Clin$msi_imputed, "CD8 T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp16 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`CD8 T cells`,
    exmat_MCP_Clin$stage, "CD8 T cells II vs III")

# En un grid 1x3 (patchwork)
bp3 + bp15 + bp16 + plot_layout(ncol = 3, nrow = 1)

5.2.3 Cytotoxic lymphocytes

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp17 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
    exmat_MCP_Clin$msi_imputed, "Cytotoxic lymphocytes MSI vs MSS")
# Diagrama de Barras II vs III
bp18 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
    exmat_MCP_Clin$stage, "Cytotoxic lymphocytes II vs III")

# En un grid 1x3 (patchwork)
bp4 + bp17 + bp18 + plot_layout(ncol = 3, nrow = 1)

5.2.4 B lineage

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp19 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`B lineage`,
    exmat_MCP_Clin$msi_imputed, "B lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp20 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`B lineage`,
    exmat_MCP_Clin$stage, "B lineage II vs III")

# En un grid 1x3 (patchwork)
bp5 + bp19 + bp20 + plot_layout(ncol = 3, nrow = 1)

5.2.5 Monocytic lineage

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp21 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Monocytic lineage`,
    exmat_MCP_Clin$msi_imputed, "Monocytic lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp22 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Monocytic lineage`,
    exmat_MCP_Clin$stage, "Monocytic lineage II vs III")

# En un grid 1x3 (patchwork)
bp6 + bp21 + bp22 + plot_layout(ncol = 3, nrow = 1)

5.2.6 Myeloid dendritic cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp23 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Myeloid dendritic cells`,
    exmat_MCP_Clin$msi_imputed, "Myeloid dendritic cells MSI vs MSS")
# Diagrama de Barras II vs III
bp24 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Myeloid dendritic cells`,
    exmat_MCP_Clin$stage, "Myeloid dendritic cells II vs III")

# En un grid 1x3 (patchwork)
bp7 + bp23 + bp24 + plot_layout(ncol = 3, nrow = 1)

5.2.7 Neutrophils

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp25 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Neutrophils,
    exmat_MCP_Clin$msi_imputed, "Neutrophils MSI vs MSS")
# Diagrama de Barras II vs III
bp26 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Neutrophils,
    exmat_MCP_Clin$stage, "Neutrophils II vs III")

# En un grid 1x3 (patchwork)
bp8 + bp25 + bp26 + plot_layout(ncol = 3, nrow = 1)

5.2.8 Endothelial cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp27 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Endothelial cells`,
    exmat_MCP_Clin$msi_imputed, "Endothelial cells MSI vs MSS")
# Diagrama de Barras II vs III
bp28 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Endothelial cells`,
    exmat_MCP_Clin$stage, "Endothelial cells II vs III")

# En un grid 1x3 (patchwork)
bp9 + bp27 + bp28 + plot_layout(ncol = 3, nrow = 1)

5.2.9 Fibroblasts

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp29 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Fibroblasts,
    exmat_MCP_Clin$msi_imputed, "Fibroblasts MSI vs MSS")
# Diagrama de Barras II vs III
bp30 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Fibroblasts,
    exmat_MCP_Clin$stage, "Fibroblasts II vs III")

# En un grid 1x3 (patchwork)
bp10 + bp29 + bp30 + plot_layout(ncol = 3, nrow = 1)

5.2.10 NK cells

# En un grid 1x3
bp1 + bp11 + bp12 + plot_layout(ncol = 3, nrow = 1)

6 Análisis CMS

La clasificación CMS divide los tumores de CRC en cuatro subtipos principales, basados en características moleculares distintas, lo que refleja diferencias en la biología del tumor, la interacción con el TME, el pronóstico y la respuesta a tratamientos específicos. (Roelands et al. 2017)

Subtipo	Características	Inmunidad	Prevalencia
CMS1 (Subtipo Inmuno)	Alta inmunogenicidad, frecuentemente asociada con MSI y una alta carga mutacional.	Tienen una infiltración activa de células inmunitarias efectoras, lo que puede hacerlos más susceptibles a terapias inmunogénicas.	Aproximadamente el 14% de los casos de CRC.
CMS2 (Subtipo Canónico)	Tumores con alta expresión de genes relacionados con el ciclo celular y la señalización WNT/MYC, presentando una marcada proliferación celular.	Son generalmente poco inmunogénicos, con menor infiltración inmune comparado con CMS1 y CMS4.	Es el subtipo más común, abarcando alrededor del 37% de los casos de CRC.
CMS3 (Subtipo Metabólico)	Presenta alteraciones metabólicas, incluyendo un metabolismo energético distintivo, con una menor carga mutacional en comparación con CMS1.	Similar a CMS2, muestra baja inmunogenicidad.	Constituye aproximadamente el 13% de los casos de CRC
CMS4 (Subtipo Mesenquimal)	Se caracteriza por la activación de vías relacionadas con la angiogénesis, la activación de células estromales y la señalización TGF-β, lo que conduce a un entorno inmunosupresor.	Aunque presenta una significativa infiltración inmune, esta es predominantemente de naturaleza supresora, con una alta presencia de células estromales.	Aproximadamente el 23% de los casos de CRC caen en este subtipo.

6.1 Infiltración CMS

| Vemos la infiltración celular según CMS e inestabilidad de microsatélites.

# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024

# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw_cms <- lapply(tipos_celulares, function(cell_type) {
    kruskal.test(reformulate("cms", response = cell_type), data = exmat_MCP_Clin)
})

# Añadimos los nombres de los tipos celulares
names(resultados_kw_cms) <- tipos_celulares

# Extraemos los valores de p de todos los tests realizados
p_valores_cms <- sapply(resultados_kw_cms, function(x) x$p.value)

# Quitamos los NA's de la columna cms
exmat_MCP_Clin_CMS <- exmat_MCP_Clin[!is.na(exmat_MCP_Clin$cms), ]
# Generamos todos los boxplot 02/04/2024
cms <- exmat_MCP_Clin_CMS$cms
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`NK cells`,
    cms, "NK cells", p_valores_cms["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`T cells`,
    cms, "T cells", p_valores_cms["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`CD8 T cells`,
    cms, "CD8 T cells", p_valores_cms["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Cytotoxic lymphocytes`,
    cms, "Cytotoxic lymphocytes", p_valores_cms["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`B lineage`,
    cms, "B lineage", p_valores_cms["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Monocytic lineage`,
    cms, "Monocytic lineage", p_valores_cms["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Myeloid dendritic cells`,
    cms, "Myeloid dendritic cells", p_valores_cms["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Neutrophils,
    cms, "Neutrophils", p_valores_cms["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Endothelial cells`,
    cms, "Endothelial cells", p_valores_cms["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Fibroblasts,
    cms, "Fibroblasts", p_valores_cms["Fibroblasts"])

# En un grid 3x2 (patchwork)
p2 + bp1 + bp2 + plot_layout(ncol = 3, nrow = 1)

bp3 + bp4 + bp5 + plot_layout(ncol = 3, nrow = 1)

bp6 + bp7 + bp8 + plot_layout(ncol = 3, nrow = 1)

bp9 + bp10 + plot_layout(ncol = 3, nrow = 1)

Los tumores CMS2 y CMS3 son poco inmunogénicos; CMS1 y CMS4 difieren en el tipo de infiltración inmune.
Los tumores CMS1 tienden a acumular una gran cantidad de mutaciones debido al estado de MSI y atraen células efectoras inmunes.
Los tumores CMS4 exhiben un TME inmunosupresor con infiltración estromal. (Lanuza et al. 2022)

6.2 Infiltración CMS-ESTADIO

| Vemos la infiltración celular según CMS y estadio.

# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 22/04/2024

# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw_cms_stage <- lapply(tipos_celulares, function(cell_type) {
    kruskal.test(reformulate("cms_stage", response = cell_type), data = exmat_MCP_Clin)
})

# Añadimos los nombres de los tipos celulares
names(resultados_kw_cms_stage) <- tipos_celulares

# Extraemos los valores de p de todos los tests realizados
p_valores_cms_stage <- sapply(resultados_kw_cms_stage, function(x) x$p.value)

# Quitamos los NA's de la columna cms
exmat_MCP_Clin_CMS_stage <- exmat_MCP_Clin[!is.na(exmat_MCP_Clin$cms_stage), ]
# Generamos todos los boxplot 02/04/2024
cms_stage <- exmat_MCP_Clin_CMS$cms_stage
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`NK cells`,
    cms_stage, "NK cells", p_valores_cms_stage["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`T cells`,
    cms_stage, "T cells", p_valores_cms_stage["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms_stage, exmat_MCP_Clin_CMS$`CD8 T cells`,
    cms_stage, "CD8 T cells", p_valores_cms_stage["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Cytotoxic lymphocytes`,
    cms_stage, "Cytotoxic lymphocytes", p_valores_cms_stage["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`B lineage`,
    cms_stage, "B lineage", p_valores_cms_stage["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Monocytic lineage`,
    cms_stage, "Monocytic lineage", p_valores_cms_stage["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Myeloid dendritic cells`,
    cms_stage, "Myeloid dendritic cells", p_valores_cms_stage["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$Neutrophils,
    cms_stage, "Neutrophils", p_valores_cms_stage["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Endothelial cells`,
    cms_stage, "Endothelial cells", p_valores_cms_stage["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$Fibroblasts,
    cms_stage, "Fibroblasts", p_valores_cms_stage["Fibroblasts"])

# En un grid 2x1 (patchwork)
bp1 + bp2 + plot_layout(ncol = 2, nrow = 1)

bp3 + bp4 + plot_layout(ncol = 2, nrow = 1)

bp5 + bp6 + plot_layout(ncol = 2, nrow = 1)

bp7 + bp8 + plot_layout(ncol = 2, nrow = 1)

bp9 + bp10 + plot_layout(ncol = 2, nrow = 1)

p5 + plot_layout(ncol = 2, nrow = 1)

7 Análisis HLA’s

Hacemos un data frame con los datos de expresión de los genes del Complejo Mayor de Hitocompatibilidad I HLA-A, HLA-B, HLA-C y el gen marcador de las células NK, NKp46 (NCR1 en nomenclatura SYMBOL). Además añadimos la columna HLA-ABC con la suma de expresión de los 3 genes.

# Nos quedamos solo con los genes del complejo mayor de hitocompatibilidad I
# HLA-ABC y con el marcador de células NK, el gen NKp46 (NCR1 en nomenclatura
# SYMBOL)
hla_nk_genes <- c("HLA-A", "HLA-B", "HLA-C", "HLA-E", "NCR1")
df_hla_nk_expression <- df_exmat_t[hla_nk_genes]
# Calculamos la mediana de la columna HLA-A
median_hla_a <- median(df_hla_nk_expression$`HLA-A`, na.rm = TRUE)
# Calculamos la mediana de la columna HLA-A
median_hla_b <- median(df_hla_nk_expression$`HLA-B`, na.rm = TRUE)
# Calculamos la mediana de la columna HLA-A
median_hla_c <- median(df_hla_nk_expression$`HLA-C`, na.rm = TRUE)
# Calculamos la mediana de la columna HLA-E
median_hla_e <- median(df_hla_nk_expression$`HLA-E`, na.rm = TRUE)

# Añadimos una nueva columna que indique si el nivel de HLA-A es 'Low' o 'High'
df_hla_nk_expression$HLA_A_cat <- ifelse(df_hla_nk_expression$`HLA-A` < median_hla_a,
    "Low", "High")
# Convertimos la columna 'HLA_A_cat' a factor
df_hla_nk_expression$HLA_A_cat <- factor(df_hla_nk_expression$HLA_A_cat, levels = c("Low",
    "High"))

# Añadimos una nueva columna que indique si el nivel de HLA-B es 'Low' o 'High'
df_hla_nk_expression$HLA_B_cat <- ifelse(df_hla_nk_expression$`HLA-B` < median_hla_b,
    "Low", "High")
# Convertimos la columna 'HLA_B_cat' a factor
df_hla_nk_expression$HLA_B_cat <- factor(df_hla_nk_expression$HLA_B_cat, levels = c("Low",
    "High"))

# Añadimos una nueva columna que indique si el nivel de HLA-C es 'Low' o 'High'
df_hla_nk_expression$HLA_C_cat <- ifelse(df_hla_nk_expression$`HLA-C` < median_hla_c,
    "Low", "High")
# Convertimos la columna 'HLA_C_cat' a factor
df_hla_nk_expression$HLA_C_cat <- factor(df_hla_nk_expression$HLA_C_cat, levels = c("Low",
    "High"))

# Añadimos una nueva columna que indique si el nivel de HLA-E es 'Low' o 'High'
df_hla_nk_expression$HLA_E_cat <- ifelse(df_hla_nk_expression$`HLA-E` < median_hla_e,
    "Low", "High")
# Convertimos la columna 'HLA_E_cat' a factor
df_hla_nk_expression$HLA_E_cat <- factor(df_hla_nk_expression$HLA_E_cat, levels = c("Low",
    "High"))

# Añadimos los datos clínicos
df_hla_nk_expression <- cbind(exmat_MCP_Clin, df_hla_nk_expression)

7.1 HLA-A

7.1.1 Boxplot HLA-A Cat vs NKp46

Boxplot comparando la expresión del gen NKp46 entre las categorías baja (Low) y alta (High) de HLA-A. Para realizar las categorías se ha utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no paramétrico de Wilcoxon para evaluar las diferencias.

# Función para generar gráficos de caja y pruebas de Mann-Whitney U
boxplot_MWU <- function(data, stage, msi, hla_var, nk_var, title_prefix) {
    # Filtrar los datos
    df_filtered <- sqldf(paste0("SELECT * FROM df_hla_nk_expression WHERE stage = '",
        stage, "' AND msi_imputed = '", msi, "'"))

    # Realizar el test de Mann-Whitney U para los grupos Low y High
    nLow <- df_filtered$NCR1[df_filtered[[hla_var]] == "Low"]
    nHigh <- df_filtered$NCR1[df_filtered[[hla_var]] == "High"]
    mwU_test <- wilcox.test(nLow, nHigh, alternative = "two.sided")

    # Generar el título completo
    title <- paste0(title_prefix, ", ", stage, ", ", msi)

    # Generar el boxplot
    bp <- boxplot_MCPcounter(data = df_filtered, x = df_filtered[[hla_var]], y = df_filtered[[nk_var]],
        fill = df_filtered[[hla_var]], title = title, p_value = mwU_test$p.value,
        xlabel = hla_var, ylabel = nk_var)

    return(bp)
}

# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp35 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_A_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp36 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_A_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp37 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_A_cat",
    nk_var = "NCR1", title_prefix = "StageIII")
bp38 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_A_cat",
    nk_var = "NCR1", title_prefix = "StageIII")

# Visualizar los gráficos en un grid 2x2 (patchwork)
bp35 + bp36 + bp37 + bp38 + plot_layout(ncol = 2, nrow = 2)

No da significación en ningún caso. No hay una diferencia signficativa en la mediana de la expresión del gen NCR1 entre los grupos con baja y alta expresión de HLA-A.

7.1.2 Estadio II - MSI/MSS HLA-A

# Seleccionamos estadio tipo II
df_hla_nk_expression_II <- sqldf('SELECT * FROM df_hla_nk_expression WHERE stage = "II"')

# Creamos el gráfico de dispersión para Estadio II
pd2 <- ggplot(df_hla_nk_expression_II, aes(x = `HLA-A`, y = NCR1)) +
  geom_point(aes(color = msi_imputed, shape = msi_imputed), size = 3, stroke = 1) +
  labs(title = "Estadio II",
       x = "HLA-A",
       y = "NKP46 (NCR1)",
       color = "MSI/MSS",
       shape = "MSI/MSS") +
  scale_shape_manual(values = c(16, 2)) +  # cuadrado lleno, triángulo
  scale_color_manual(values = c("darkorange2", "dodgerblue4")) +
  theme_classic() +
  theme(
    axis.line = element_line(color = "black"),
    legend.position = c(0.1, 0.9))

# Añadimos líneas de regresión lineal para cada grupo MSI y MSS
pd2 <- pd2 + geom_smooth(method = "lm", aes(color = msi_imputed), se = FALSE)

# Visualizamos el gráfico
pd2

# Ajustamos un modelo lineal para Estadio II
lmod_ser <- lm(NCR1 ~ `HLA-A` * msi_imputed, data = df_hla_nk_expression_II)
# Hacemos una contraste de paralelismo de rectas con anova
anova(lmod_ser)

Analysis of Variance Table

Response: NCR1
                     Df Sum Sq Mean Sq F value   Pr(>F)   
`HLA-A`               1   0.16  0.1625  0.2164 0.641981   
msi_imputed           1   1.49  1.4871  1.9799 0.159859   
`HLA-A`:msi_imputed   1   6.93  6.9312  9.2282 0.002474 **
Residuals           680 510.74  0.7511                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El efecto principal de HLA-A sobre NCR1 no es significativo (p = 0.641981). Esto sugiere que no hay suficiente evidencia para concluir que HLA-A tiene un efecto significativo en NCR1 por sí solo.

El efecto msi_imputed sobre NCR1 no es significativo (p = 0.159859). Esto indica que no hay suficiente evidencia para concluir que msi_imputed tiene un efecto significativo en NCR1 por sí solo.

Vemos que hay significación en la hipótesis de paralelismo (Pr: 0.002474).La interacción entre HLA-A y msi_imputed es significativa. Esto indica que el efecto de HLA-A sobre NCR1 depende de los valores de msi_imputed. Dado que el p-valor es menor que 0.01 (**), hay una fuerte evidencia contra la hipótesis nula de que no hay interacción entre HLA-A y msi_imputed.

Aquí se observa para los tumores MSI una relación inversa entre la expresión de NCR1 y HLA-A; Cuanta menor expresión de HLA’s hay una mayor expresión de NKp46, es decir, cuanta menor expresión de HLA-A hay una mayor infiltración de células NK activadas.

7.1.3 Estadio III - MSI/MSS HLA-A

# Seleccionamos estadio tipo II
df_hla_nk_expression_III <- sqldf('SELECT * FROM df_hla_nk_expression WHERE stage = "III"')

# Creamos el gráfico de dispersión para Estadio II
pd3 <- ggplot(df_hla_nk_expression_III, aes(x = `HLA-A`, y = NCR1)) +
  geom_point(aes(color = msi_imputed, shape = msi_imputed), size = 3, stroke = 1) +
  labs(title = "Estadio III",
       x = "HLA-A",
       y = "NKP46 (NCR1)",
       color = "MSI/MSS",
       shape = "MSI/MSS") +
  scale_shape_manual(values = c(16, 2)) +  # cuadrado lleno, triángulo
  scale_color_manual(values = c("darkorange2", "dodgerblue4")) +
  theme_classic() +
  theme(
    axis.line = element_line(color = "black"),
    legend.position = c(0.1, 0.9))

# Añadimos líneas de regresión lineal para cada grupo MSI y MSS
pd3 <- pd3 + geom_smooth(method = "lm", aes(color = msi_imputed), se = FALSE)

# Visualizamos el gráfico
pd3

# Ajustamos un modelo lineal para Estadio II
lmod_ser <- lm(NCR1 ~ `HLA-A` * msi_imputed, data = df_hla_nk_expression_III)
# Hacemos una contraste de paralelismo de rectas con anova
anova(lmod_ser)

Analysis of Variance Table

Response: NCR1
                     Df  Sum Sq Mean Sq F value   Pr(>F)   
`HLA-A`               1   3.172  3.1717  4.8625 0.028028 * 
msi_imputed           1   5.148  5.1481  7.8925 0.005213 **
`HLA-A`:msi_imputed   1   0.036  0.0356  0.0545 0.815465   
Residuals           391 255.037  0.6523                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El efecto principal de HLA-A sobre NCR1 es significativo (p = 0.028028). Esto sugiere que hay una relación significativa entre HLA-A y NCR1 por sí solo. Dado que el p-valor es menor que 0.05 (*), hay suficiente evidencia para rechazar la hipótesis nula de que HLA-A no tiene un efecto sobre NCR1.

El efecto principal de msi_imputed sobre NCR1 es significativo (p = 0.005213). Esto indica que hay una relación significativa entre msi_imputed y NCR1 por sí solo. Dado que el p-valor es menor que 0.01 (**), hay una fuerte evidencia para rechazar la hipótesis nula de que msi_imputed no tiene un efecto sobre NCR1.

La interacción entre HLA-A y msi_imputed no es significativa (p = 0.815465). Esto indica que no hay suficiente evidencia para concluir que el efecto de HLA-A sobre NCR1 depende de los valores de msi_imputed.

En este caso no hay significación en la hipótesis de paralelismo (Pr: 0.738591) aunque sí en la de coincidencia de las rectas (Pr: 0.006508). Resultado también de difícil interpretación: los efectos principales de HLA-ABC y msi_imputed son ambos estadísticamente significativos y aportan individualmente a la explicación de la variabilidad en NCR1, con msi_imputed mostrando un efecto algo más fuerte según valor p más bajos.

Aquí se intuye una relación inversa entre la expresión de NCR1 y HLA-A; es decir, una relación inversa entre infiltración de células NK y expresión de HLA’s.

7.2 HLA-B

7.2.1 Boxplot HLA-B Cat vs NKp46

Boxplot comparando la expresión del gen NKp46 entre las categorías baja (Low) y alta (High) de HLA-B. Para realizar las categorías se ha utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no paramétrico de Wilcoxon para evaluar las diferencias.

# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp39 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_B_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp40 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_B_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp41 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_B_cat",
    nk_var = "NCR1", title_prefix = "StageIII")
bp42 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_B_cat",
    nk_var = "NCR1", title_prefix = "StageIII")


# Visualizar los gráficos en un grid 2x2 (patchwork)
bp39 + bp40 + bp41 + bp42 + plot_layout(ncol = 2, nrow = 2)

Nos da significación en el primer caso StageII, MSI.

7.3 HLA-C

7.3.1 Boxplot HLA-C Cat vs NKp46

Boxplot comparando la expresión del gen NKp46 entre las categorías baja (Low) y alta (High) de HLA-C. Para realizar las categorías se ha utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no paramétrico de Wilcoxon para evaluar las diferencias.

# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp43 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_C_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp44 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_C_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp45 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_C_cat",
    nk_var = "NCR1", title_prefix = "StageIII")
bp46 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_C_cat",
    nk_var = "NCR1", title_prefix = "StageIII")

# Visualizar los gráficos en un grid 2x2 (patchwork)
bp43 + bp44 + bp45 + bp46 + plot_layout(ncol = 2, nrow = 2)

Nos da significación en los casos StageII, MSS y StageIII, MSS.

7.4 HLA-E

7.4.1 Boxplot HLA-E Cat vs NKp46

Boxplot comparando la expresión del gen NKp46 entre las categorías baja (Low) y alta (High) de HLA-E. Para realizar las categorías se ha utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no paramétrico de Wilcoxon para evaluar las diferencias.

bp47 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_E_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp48 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_E_cat",
    nk_var = "NCR1", title_prefix = "StageII")
bp49 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_E_cat",
    nk_var = "NCR1", title_prefix = "StageIII")
bp50 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_E_cat",
    nk_var = "NCR1", title_prefix = "StageIII")

# Visualizar los gráficos en un grid 2x2
bp47 + bp48 + bp49 + bp50 + plot_layout(ncol = 2, nrow = 2)

Nos da significación en los casos StageII, MSS y StageIII, MSS.

7.4.2 MSS Estadio II/III HLA-E

# Seleccionamos estadio tipo II ya lo hemos hecho antes
df_hla_nk_expression_MSS <- sqldf('SELECT * FROM df_hla_nk_expression WHERE msi_imputed = "MSS"')

# Creamos el gráfico de dispersión para Estadio II
pd4 <- ggplot(df_hla_nk_expression_MSS, aes(x = `HLA-E`, y = NCR1)) +
  geom_point(aes(color = stage, shape = stage), size = 3, stroke = 1) +
  labs(title = "MSS",
       x = "HLA-E",
       y = "NKP46 (NCR1)",
       color = "II/III",
       shape = "II/III") +
  scale_shape_manual(values = c(16, 2)) +  # cuadrado lleno, triángulo
  scale_color_manual(values = c("darkorange2", "dodgerblue4")) +
  theme_classic() +
  theme(
    axis.line = element_line(color = "black"),
    legend.position = c(0.1, 0.9))

# Añadimos líneas de regresión lineal para cada grupo MSI y MSS
pd4 <- pd4 + geom_smooth(method = "lm", aes(color = stage), se = FALSE)

# Visualizamos el gráfico
pd4

# Ajustamos un modelo lineal para Estadio II
lmod_ser <- lm(NCR1 ~ `HLA-A` * stage, data = df_hla_nk_expression_MSS)
# Hacemos una contraste de paralelismo de rectas con anova
anova(lmod_ser)

Analysis of Variance Table

Response: NCR1
               Df Sum Sq Mean Sq F value   Pr(>F)   
`HLA-A`         1   5.67  5.6717  8.1903 0.004317 **
stage           1   1.23  1.2299  1.7760 0.183005   
`HLA-A`:stage   1   0.00  0.0013  0.0018 0.965735   
Residuals     836 578.92  0.6925                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El efecto principal de HLA-A en NCR1 es significativo (p = 0.004317). Esto significa que hay una relación significativa entre HLA-A y la infiltración de células NK activadas. Dado que el p-valor es menor que 0.01 (**), esto indica una fuerte evidencia contra la hipótesis nula de que HLA-A no tiene un efecto sobre NCR1.

El efecto de stage en NCR1 no es significativo (p = 0.183005). Esto sugiere que no hay suficiente evidencia para concluir que el estadio tenga un efecto sobre NCR1 cuando se considera solo este factor.

Vemos que hay significación en la hipótesis de paralelismo HLA-A:stage (Pr: 0.965735), no son paralelas. La interacción entre HLA-A y stage no es significativa (p = 0.965735). Esto indica que no hay suficiente evidencia para concluir que el efecto de HLA-A sobre NCR1 depende del estadio de la muestra.

8 Análisi HALA’s con infiltración (MCPcounter)

8.1 Boxplot HLA-A Cat vs ‘NK cells’

# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp51 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_A_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp52 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_A_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp53 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_A_cat",
    nk_var = "NK cells", title_prefix = "StageIII")
bp54 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_A_cat",
    nk_var = "NK cells", title_prefix = "StageIII")

# Visualizar los gráficos en un grid 2x2 (patchwork)
bp51 + bp52 + bp53 + bp54 + plot_layout(ncol = 2, nrow = 2)

8.2 Boxplot HLA-B Cat vs ‘NK cells’

# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp55 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_B_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp56 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_B_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp57 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_B_cat",
    nk_var = "NK cells", title_prefix = "StageIII")
bp58 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_B_cat",
    nk_var = "NK cells", title_prefix = "StageIII")

# Visualizar los gráficos en un grid 2x2 (patchwork)
bp55 + bp56 + bp57 + bp58 + plot_layout(ncol = 2, nrow = 2)

8.3 Boxplot HLA-C Cat vs ‘NK cells’

# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp59 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_C_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp60 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_C_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp61 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_C_cat",
    nk_var = "NK cells", title_prefix = "StageIII")
bp62 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_C_cat",
    nk_var = "NK cells", title_prefix = "StageIII")

# Visualizar los gráficos en un grid 2x2 (patchwork)
bp59 + bp60 + bp61 + bp62 + plot_layout(ncol = 2, nrow = 2)

8.4 Boxplot HLA-E Cat vs ‘NK cells’

# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp63 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_E_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp64 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_E_cat",
    nk_var = "NK cells", title_prefix = "StageII")
bp65 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_E_cat",
    nk_var = "NK cells", title_prefix = "StageIII")
bp66 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_E_cat",
    nk_var = "NK cells", title_prefix = "StageIII")

# Visualizar los gráficos en un grid 2x2 (patchwork)
bp63 + bp64 + bp65 + bp66 + plot_layout(ncol = 2, nrow = 2)

9 Algoritmo Random Forest (ML)

El algoritmo “Random Forest” es un método de aprendizaje supervisado utilizado para la clasificación en el campo de machine learning (ML). Se utilizan en una amplia gama de campos como, por ejemplo, en el diagnóstico médico.
Aplicaremos este algoritmo a nuestros datos para clasificar los CRC así como el estadio de los mismos.

# library(randomForest) library(caret) library(dplyr)

# Función para crear y evaluar un modelo Random Forest con 100 árboles
run_randomForest_model <- function(data, target_name = "class", p_train = 0.67, seed = 12357,
    ntree = 100) {
    set.seed(seed)

    # Preparamos los datos
    clean_data <- na.omit(data)
    # Obtenemos el número total de filas (datos) n <- nrow(clean_data)
    # Obtenemos los índices de filas al azar para entrenamiento indices_train
    # <- sample(1:n, round(p_train * n))

    # Usar createDataPartition para obtener índices de entrenamiento
    indices_train <- createDataPartition(y = clean_data$class, p = params$p.train,
        list = FALSE)

    # Seleccionamos las filas para entrenamiento
    data_train <- clean_data[indices_train, ]
    # Seleccionamos las filas restantes para pruebas
    data_test <- clean_data[-indices_train, ]
    # Creamos los labels de class de 'train'
    train_labels <- data_train[[target_name]]
    # Creamos los labels de class de 'test'
    test_labels <- data_test[[target_name]]

    # Creamos el modelo
    model <- randomForest(data_train[-1], train_labels, ntree = ntree)

    # Predicción y evaluación
    predictions <- predict(model, data_test[-1], type = "response")
    results <- confusionMatrix(table(predictions, test_labels))

    return(list(model = model, results = results))
}

9.1 Estadio

# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_stage <- cbind(df_dataClin["stage"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_stage)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = 100) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 37.36%
Confusion matrix:
     II III class.error
II  413  43  0.09429825
III 226  38  0.85606061

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions  II III
        II  215 114
        III  13  17
                                          
               Accuracy : 0.6462          
                 95% CI : (0.5943, 0.6957)
    No Information Rate : 0.6351          
    P-Value [Acc > NIR] : 0.3522          
                                          
                  Kappa : 0.087           
                                          
 Mcnemar's Test P-Value : <2e-16          
                                          
            Sensitivity : 0.9430          
            Specificity : 0.1298          
         Pos Pred Value : 0.6535          
         Neg Pred Value : 0.5667          
             Prevalence : 0.6351          
         Detection Rate : 0.5989          
   Detection Prevalence : 0.9164          
      Balanced Accuracy : 0.5364          
                                          
       'Positive' Class : II              
                                          

# Regogemos los datos para comparar de la predicción
stage_Accuracy <- model_results$results$overall["Accuracy"]
stage_Kappa <- model_results$results$overall["Kappa"]
stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.2 MSI/MSS

# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_msi <- cbind(df_dataClin["msi_imputed"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_msi)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_msi, ntree = 200)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 200
No. of variables tried at each split: 123

        OOB estimate of  error rate: 12.64%
Confusion matrix:
    MSI MSS class.error
MSI  94  66  0.41250000
MSS  25 535  0.04464286

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions MSI MSS
        MSI  34   7
        MSS  45 273
                                          
               Accuracy : 0.8552          
                 95% CI : (0.8144, 0.8899)
    No Information Rate : 0.7799          
    P-Value [Acc > NIR] : 0.0002111       
                                          
                  Kappa : 0.49            
                                          
 Mcnemar's Test P-Value : 2.882e-07       
                                          
            Sensitivity : 0.43038         
            Specificity : 0.97500         
         Pos Pred Value : 0.82927         
         Neg Pred Value : 0.85849         
             Prevalence : 0.22006         
         Detection Rate : 0.09471         
   Detection Prevalence : 0.11421         
      Balanced Accuracy : 0.70269         
                                          
       'Positive' Class : MSI             
                                          

# Regogemos los datos para comparar de la predicción
msi_Accuracy <- model_results$results$overall["Accuracy"]
msi_Kappa <- model_results$results$overall["Kappa"]
msi_AccuracyLower <- model_results$results$overall["AccuracyLower"]
msi_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.3 MSI/MSS-Estadio

# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_msi_stage <- cbind(df_dataClin["grupo"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_msi_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_msi_stage, ntree = 200)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 200
No. of variables tried at each split: 123

        OOB estimate of  error rate: 45.08%
Confusion matrix:
        II_MSI II_MSS III_MSI III_MSS class.error
II_MSI      71     43       0       2   0.3879310
II_MSS      13    291       0      37   0.1466276
III_MSI     26     12       1       5   0.9772727
III_MSS      4    183       0      33   0.8500000

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions II_MSI II_MSS III_MSI III_MSS
    II_MSI      34      8      13       8
    II_MSS      22    154       8      91
    III_MSI      0      0       0       0
    III_MSS      1      8       1      10

Overall Statistics
                                          
               Accuracy : 0.5531          
                 95% CI : (0.4999, 0.6053)
    No Information Rate : 0.4749          
    P-Value [Acc > NIR] : 0.001813        
                                          
                  Kappa : 0.2428          
                                          
 Mcnemar's Test P-Value : < 2.2e-16       

Statistics by Class:

                     Class: II_MSI Class: II_MSS Class: III_MSI Class: III_MSS
Sensitivity                0.59649        0.9059        0.00000        0.09174
Specificity                0.90365        0.3564        1.00000        0.95984
Pos Pred Value             0.53968        0.5600            NaN        0.50000
Neg Pred Value             0.92203        0.8072        0.93855        0.70710
Prevalence                 0.15922        0.4749        0.06145        0.30447
Detection Rate             0.09497        0.4302        0.00000        0.02793
Detection Prevalence       0.17598        0.7682        0.00000        0.05587
Balanced Accuracy          0.75007        0.6311        0.50000        0.52579

# Regogemos los datos para comparar de la predicción
msi_mss_stage_Accuracy <- model_results$results$overall["Accuracy"]
msi_mss_stage_Kappa <- model_results$results$overall["Kappa"]
msi_mss_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
msi_mss_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.4 CMS

# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_cms <- cbind(df_dataClin["cms"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = 100) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 11.36%
Confusion matrix:
     CMS1 CMS2 CMS3 CMS4 class.error
CMS1  109    7    4   10  0.16153846
CMS2    0  214    2    4  0.02727273
CMS3    9   20   61    1  0.32967033
CMS4    4    9    0  162  0.07428571

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions CMS1 CMS2 CMS3 CMS4
       CMS1   58    0    1    2
       CMS2    1  106    9    6
       CMS3    2    1   35    0
       CMS4    3    2    0   79

Overall Statistics
                                          
               Accuracy : 0.9115          
                 95% CI : (0.8738, 0.9409)
    No Information Rate : 0.3574          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.8767          
                                          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: CMS1 Class: CMS2 Class: CMS3 Class: CMS4
Sensitivity               0.9062      0.9725      0.7778      0.9080
Specificity               0.9876      0.9184      0.9885      0.9771
Pos Pred Value            0.9508      0.8689      0.9211      0.9405
Neg Pred Value            0.9754      0.9836      0.9625      0.9638
Prevalence                0.2098      0.3574      0.1475      0.2852
Detection Rate            0.1902      0.3475      0.1148      0.2590
Detection Prevalence      0.2000      0.4000      0.1246      0.2754
Balanced Accuracy         0.9469      0.9454      0.8831      0.9426

# Regogemos los datos para comparar de la predicción
cms_Accuracy <- model_results$results$overall["Accuracy"]
cms_Kappa <- model_results$results$overall["Kappa"]
cms_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.5 CMS-Estadio

# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_cms_stage <- cbind(df_dataClin["cms_stage"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms_stage, ntree = 200)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 200
No. of variables tried at each split: 123

        OOB estimate of  error rate: 44.16%
Confusion matrix:
         II_CMS1 II_CMS2 II_CMS3 II_CMS4 III_CMS1 III_CMS2 III_CMS3 III_CMS4
II_CMS1       66       1       4       8        3        1        0        0
II_CMS2        0     124       1       4        0       11        0        0
II_CMS3        9      11      43       0        0        1        0        0
II_CMS4        0       7       0      78        0        0        0       16
III_CMS1      28       5       1       5        7        0        0        1
III_CMS2       0      71       2       1        0        6        0        0
III_CMS3       1       6      16       2        0        0        2        0
III_CMS4       0       7       0      47        0        2        0       18
         class.error
II_CMS1    0.2048193
II_CMS2    0.1142857
II_CMS3    0.3281250
II_CMS4    0.2277228
III_CMS1   0.8510638
III_CMS2   0.9250000
III_CMS3   0.9259259
III_CMS4   0.7567568

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions II_CMS1 II_CMS2 II_CMS3 II_CMS4 III_CMS1 III_CMS2 III_CMS3 III_CMS4
   II_CMS1       38       0       2       1       17        0        0        0
   II_CMS2        0      66       8       2        2       35        2        1
   II_CMS3        0       0      22       0        1        0        9        0
   II_CMS4        2       2       0      42        2        0        0       30
   III_CMS1       1       0       0       0        1        0        0        0
   III_CMS2       0       2       0       0        0        4        2        0
   III_CMS3       0       0       0       0        0        0        0        0
   III_CMS4       0       0       0       5        0        0        0        6

Overall Statistics
                                          
               Accuracy : 0.5869          
                 95% CI : (0.5294, 0.6427)
    No Information Rate : 0.2295          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.4999          
                                          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: II_CMS1 Class: II_CMS2 Class: II_CMS3
Sensitivity                  0.9268         0.9429        0.68750
Specificity                  0.9242         0.7872        0.96337
Pos Pred Value               0.6552         0.5690        0.68750
Neg Pred Value               0.9879         0.9788        0.96337
Prevalence                   0.1344         0.2295        0.10492
Detection Rate               0.1246         0.2164        0.07213
Detection Prevalence         0.1902         0.3803        0.10492
Balanced Accuracy            0.9255         0.8650        0.82543
                     Class: II_CMS4 Class: III_CMS1 Class: III_CMS2
Sensitivity                  0.8400        0.043478         0.10256
Specificity                  0.8588        0.996454         0.98496
Pos Pred Value               0.5385        0.500000         0.50000
Neg Pred Value               0.9648        0.927393         0.88215
Prevalence                   0.1639        0.075410         0.12787
Detection Rate               0.1377        0.003279         0.01311
Detection Prevalence         0.2557        0.006557         0.02623
Balanced Accuracy            0.8494        0.519966         0.54376
                     Class: III_CMS3 Class: III_CMS4
Sensitivity                  0.00000         0.16216
Specificity                  1.00000         0.98134
Pos Pred Value                   NaN         0.54545
Neg Pred Value               0.95738         0.89456
Prevalence                   0.04262         0.12131
Detection Rate               0.00000         0.01967
Detection Prevalence         0.00000         0.03607
Balanced Accuracy            0.50000         0.57175

# Regogemos los datos para comparar de la predicción
cms_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms_stage_Kappa <- model_results$results$overall["Kappa"]
cms_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.6 MSI-Estadio

# Filtramos y escogemos solo los datos msi_imputed = MSI
df_dataClin_exmat_msi <- df_dataClin_exmat_t[df_dataClin_exmat_t$msi_imputed == "MSI",
    ]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_msi_stage <- df_dataClin_exmat_msi %>%
    select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_msi_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_msi_stage)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 32.5%
Confusion matrix:
     II III class.error
II  106  10   0.0862069
III  42   2   0.9545455

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions II III
        II  57  22
        III  0   0
                                          
               Accuracy : 0.7215          
                 95% CI : (0.6093, 0.8165)
    No Information Rate : 0.7215          
    P-Value [Acc > NIR] : 0.5571          
                                          
                  Kappa : 0               
                                          
 Mcnemar's Test P-Value : 7.562e-06       
                                          
            Sensitivity : 1.0000          
            Specificity : 0.0000          
         Pos Pred Value : 0.7215          
         Neg Pred Value :    NaN          
             Prevalence : 0.7215          
         Detection Rate : 0.7215          
   Detection Prevalence : 1.0000          
      Balanced Accuracy : 0.5000          
                                          
       'Positive' Class : II              
                                          

# Regogemos los datos para comparar de la predicción
msi_stage_Accuracy <- model_results$results$overall["Accuracy"]
msi_stage_Kappa <- model_results$results$overall["Kappa"]
msi_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
msi_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.7 MSS-Estadio

# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_mss <- df_dataClin_exmat_t[df_dataClin_exmat_t$msi_imputed == "MSS",
    ]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_mss_stage <- df_dataClin_exmat_mss %>%
    select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_mss_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_mss_stage)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 39.57%
Confusion matrix:
     II III class.error
II  287  54   0.1583578
III 168  52   0.7636364

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions  II III
        II  155  86
        III  15  23
                                          
               Accuracy : 0.638           
                 95% CI : (0.5786, 0.6944)
    No Information Rate : 0.6093          
    P-Value [Acc > NIR] : 0.1789          
                                          
                  Kappa : 0.139           
                                          
 Mcnemar's Test P-Value : 3.278e-12       
                                          
            Sensitivity : 0.9118          
            Specificity : 0.2110          
         Pos Pred Value : 0.6432          
         Neg Pred Value : 0.6053          
             Prevalence : 0.6093          
         Detection Rate : 0.5556          
   Detection Prevalence : 0.8638          
      Balanced Accuracy : 0.5614          
                                          
       'Positive' Class : II              
                                          

# Regogemos los datos para comparar de la predicción
mss_stage_Accuracy <- model_results$results$overall["Accuracy"]
mss_stage_Kappa <- model_results$results$overall["Kappa"]
mss_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
mss_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.8 CMS1-Estadio

# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms1 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS1",
    ]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms1_stage <- df_dataClin_exmat_cms1 %>%
    select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms1_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms1_stage)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 36.15%
Confusion matrix:
    II III class.error
II  73  10   0.1204819
III 37  10   0.7872340

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions II III
        II  39  19
        III  2   4
                                          
               Accuracy : 0.6719          
                 95% CI : (0.5431, 0.7841)
    No Information Rate : 0.6406          
    P-Value [Acc > NIR] : 0.3520305       
                                          
                  Kappa : 0.1494          
                                          
 Mcnemar's Test P-Value : 0.0004803       
                                          
            Sensitivity : 0.9512          
            Specificity : 0.1739          
         Pos Pred Value : 0.6724          
         Neg Pred Value : 0.6667          
             Prevalence : 0.6406          
         Detection Rate : 0.6094          
   Detection Prevalence : 0.9062          
      Balanced Accuracy : 0.5626          
                                          
       'Positive' Class : II              
                                          

# Regogemos los datos para comparar de la predicción
cms1_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms1_stage_Kappa <- model_results$results$overall["Kappa"]
cms1_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms1_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.9 CMS2-Estadio

# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms2 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS2",
    ]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms2_stage <- df_dataClin_exmat_cms2 %>%
    select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms2_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms2_stage)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 40.45%
Confusion matrix:
     II III class.error
II  124  16   0.1142857
III  73   7   0.9125000

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions II III
        II  65  36
        III  5   3
                                         
               Accuracy : 0.6239         
                 95% CI : (0.526, 0.7148)
    No Information Rate : 0.6422         
    P-Value [Acc > NIR] : 0.6937         
                                         
                  Kappa : 0.0067         
                                         
 Mcnemar's Test P-Value : 2.797e-06      
                                         
            Sensitivity : 0.92857        
            Specificity : 0.07692        
         Pos Pred Value : 0.64356        
         Neg Pred Value : 0.37500        
             Prevalence : 0.64220        
         Detection Rate : 0.59633        
   Detection Prevalence : 0.92661        
      Balanced Accuracy : 0.50275        
                                         
       'Positive' Class : II             
                                         

# Regogemos los datos para comparar de la predicción
cms2_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms2_stage_Kappa <- model_results$results$overall["Kappa"]
cms2_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms2_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.10 CMS3-Estadio

# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms3 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS3",
    ]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms3_stage <- df_dataClin_exmat_cms3 %>%
    select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms3_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms3_stage)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 28.57%
Confusion matrix:
    II III class.error
II  63   1   0.0156250
III 25   2   0.9259259

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions II III
        II  31  13
        III  1   0
                                          
               Accuracy : 0.6889          
                 95% CI : (0.5335, 0.8183)
    No Information Rate : 0.7111          
    P-Value [Acc > NIR] : 0.695081        
                                          
                  Kappa : -0.043          
                                          
 Mcnemar's Test P-Value : 0.003283        
                                          
            Sensitivity : 0.9688          
            Specificity : 0.0000          
         Pos Pred Value : 0.7045          
         Neg Pred Value : 0.0000          
             Prevalence : 0.7111          
         Detection Rate : 0.6889          
   Detection Prevalence : 0.9778          
      Balanced Accuracy : 0.4844          
                                          
       'Positive' Class : II              
                                          

# Regogemos los datos para comparar de la predicción
cms3_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms3_stage_Kappa <- model_results$results$overall["Kappa"]
cms3_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms3_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.11 CMS4-Estadio

# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms4 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS4",
    ]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms4_stage <- df_dataClin_exmat_cms4 %>%
    select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms4_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms4_stage)

# Vemos la información del modelo
model_results$model

Call:
 randomForest(x = data_train[-1], y = train_labels, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 123

        OOB estimate of  error rate: 42.86%
Confusion matrix:
    II III class.error
II  76  25   0.2475248
III 50  24   0.6756757

# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)

Confusion Matrix and Statistics

           test_labels
predictions II III
        II  42  30
        III  8   7
                                          
               Accuracy : 0.5632          
                 95% CI : (0.4526, 0.6694)
    No Information Rate : 0.5747          
    P-Value [Acc > NIR] : 0.6292447       
                                          
                  Kappa : 0.0316          
                                          
 Mcnemar's Test P-Value : 0.0006577       
                                          
            Sensitivity : 0.8400          
            Specificity : 0.1892          
         Pos Pred Value : 0.5833          
         Neg Pred Value : 0.4667          
             Prevalence : 0.5747          
         Detection Rate : 0.4828          
   Detection Prevalence : 0.8276          
      Balanced Accuracy : 0.5146          
                                          
       'Positive' Class : II              
                                          

# Regogemos los datos para comparar de la predicción
cms4_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms4_stage_Kappa <- model_results$results$overall["Kappa"]
cms4_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms4_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]

9.12 Resultados

Modelo	Accuraccy	kappa	95% CI
Estadio	0.6462396	0.0870262	(0.5943389 - 0.6957054)
MSI/MSS	0.8551532	0.4899732	(0.8144282 - 0.8898921)
MSI/MSS-Estadio	0.5530726	0.2427587	(0.4999164 - 0.6053439)
CMS	0.9114754	0.8767308	(0.8738159 - 0.9408505)
CMS-Estadio	0.5868852	0.4998568	(0.5293643 - 0.6427016)
MSI-Estadio	0.721519	0	(0.6092649 - 0.816545)
MSS-Estadio	0.6379928	0.1390204	(0.5785781 - 0.6944416)
CMS1-Estadio	0.671875	0.1493671	(0.5431232 - 0.784128)
CMS2-Estadio	0.6238532	0.0066681	(0.525959 - 0.7148372)
CMS3-Estadio	0.6888889	-0.0430464	(0.533509 - 0.8183412)
CMS4-Estadio	0.5632184	0.0316344	(0.4526478 - 0.669361)

10 Referencias

Repositorio: (Ferrón 2024)

Becht, Etienne, Nicolas A Giraldo, Laetitia Lacroix, Bénédicte Buttard, Nabila Elarouci, Florent Petitprez, Janick Selves, et al. 2016. “Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression.” Genome Biology 17: 1–20.

Ferrón, Sheddad Kaid-Salah. 2024. “Caracterización Inmunonómica Del Microambiente Tumoral En Cáncer Colorrectal En Estadios Tempranos y Tardı́os.” https://github.com/Sheddad/TFM_Inmunonomica_TME_CRC.

Jorissen, Robert N, Peter Gibbs, Michael Christie, Saurabh Prakash, Lara Lipton, Jayesh Desai, David Kerr, et al. 2009. “Metastasis-Associated Gene Expression Changes Predict Poor Outcomes in Patients with Dukes Stage b and c Colorectal Cancer.” Clinical Cancer Research 15 (24): 7642–51.

Jorissen, Robert N, Lara Lipton, Peter Gibbs, Matthew Chapman, Jayesh Desai, Ian T Jones, Timothy J Yeatman, et al. 2008. “DNA Copy-Number Alterations Underlie Gene Expression Differences Between Microsatellite Stable and Unstable Colorectal Cancers.” Clinical Cancer Research 14 (24): 8061–69.

Lanuza, Pilar M, M Henar Alonso, Iratxe Uranga-Murillo, Sandra Garcı́a-Mulero, Cecilia Pesini, Eva M Gálvez, Ariel Ramı́rez-Labrada, Rebeca Sanz-Pamplona, and Julián Pardo. 2022. “Adoptive NK Cell Transfer as a Treatment in Colorectal Cancer Patients: Analyses of Tumour Cell Determinants Correlating with Efficacy in Vitro and in Vivo.” Frontiers in Immunology 13: 890836.

Marisa, Laetitia, Aurélien de Reyniès, Alex Duval, Janick Selves, Marie Pierre Gaub, Laure Vescovo, Marie-Christine Etienne-Grimaldi, et al. 2013. “Gene Expression Classification of Colon Cancer into Molecular Subtypes: Characterization, Validation, and Prognostic Value.” PLoS Medicine 10 (5): e1001453.

Roelands, Jessica, Peter JK Kuppen, Louis Vermeulen, Cristina Maccalli, Julie Decock, Ena Wang, Francesco M Marincola, Davide Bedognetti, and Wouter Hendrickx. 2017. “Immunogenomic Classification of Colorectal Cancer and Therapeutic Implications.” International Journal of Molecular Sciences 18 (10): 2229.

Smith, J Joshua, Natasha G Deane, Fei Wu, Nipun B Merchant, Bing Zhang, Aixiang Jiang, Pengcheng Lu, et al. 2010. “Experimentally Derived Metastasis Gene Expression Profile Predicts Recurrence and Death in Patients with Colon Cancer.” Gastroenterology 138 (3): 958–68.